home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / fitls / fitls.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  8.9 KB  |  387 lines
  1. /*
  2.     Ask for 
  3.         Data file, columns    a.dat,1,2
  4.         Equation        sin(a*x)+x^2+b+c*x^3
  5.         Initial values for a,b,c,...
  6.         Output file name
  7.     Write out results, fit, gle file.
  8. */        
  9. #include "all.h"
  10. #include <math.h>
  11. #include <time.h>
  12. int load_data(char *s, int c1, int c2);
  13. #ifdef unix
  14. #define VAXC 1
  15. #endif
  16. #ifdef VAXC
  17. #ifndef unix
  18. #define unlink delete
  19. #endif
  20. #define farfree free
  21. #define farmalloc malloc
  22. #else
  23. #include <alloc.h>
  24. #endif
  25. typedef double real;
  26.  
  27. int stringlower(char *s);
  28. int build_equation(char *expr, double pms[], int naz);
  29. int stripspace(char *ss);
  30. int var_find_az(int v_idx[],int *naz);
  31. int v_eval(double *v);
  32. /*---------------------------------------------------------------------------*/
  33. char outfile[255],inputfile[255];
  34. double anstol=1e-4;
  35. char equation[255];
  36. int pnt_alloc(int size);
  37. int col1,col2;
  38. double xmin,xmax;
  39. int gle_debug=0;
  40. #define dbg if (gle_debug>0)
  41. extern char file_name[];
  42. main()
  43. {
  44.     char eqn[200],buff[200],myfile[200],lastone[200];
  45.     char lastdat[200];
  46.     int c1,c2;
  47.     char *s1,*s2;
  48.  
  49.     lastone_get(lastone,lastdat);
  50.     printf("Input data file (x and y columns optional) [%s] ? ",lastdat); 
  51.     gets(myfile);
  52.     if (strlen(myfile)==0) strcpy(myfile,lastdat);
  53.     strcpy(lastdat,myfile);
  54.     s1 = strchr(myfile,',');
  55.     c1 = 1; c2 = 2;
  56.     if (s1!=NULL) {
  57.         s2 = strchr(s1+1,',');
  58.         c1 = atoi(s1+1);
  59.         if (s2!=NULL) c2 = atoi(s2+1);
  60.         if (c1==0 || c2==0) printf("Expecting filename.dat,xcol,ycol\n");
  61.         *s1 = 0;
  62.     }
  63.     printf("Loading data from file, %s,  xcolumn=%d, ycolumn=%d\n",myfile,c1,c2);
  64.     col1 = c1; col2 = c2; strcpy(inputfile,myfile);
  65.     load_data(myfile,c1,c2);    
  66.     printf("Valid operators: +, -, *, /, ^ (power) \n");
  67.     printf("Valid functions:\n");
  68.     printf("\t abs(), atn(), cos(), exp(), fix(), int()\n");
  69.     printf("\t log(), log10(), not(), rnd(), sgn(), sin()\n");
  70.     printf("\t sqr(), sqrt(), tan()\n");
  71.     printf("\n Enter a function of 'x' using constants 'a'...'z' \n");
  72.     printf(" e.g.\t a + b*x   (standard linear least squares fit) \n");
  73.     printf("     \t sin(x)*a+b \n");
  74.     printf("     \t a + b*x + c*x^2 + d*x^3  \n");
  75.     printf("     \t log(a*x)+(b+x)*c+a \n");
  76.     printf("\nEquation [%s]? ",lastone);
  77.     gets(eqn);
  78.     if (strlen(eqn)==0) strcpy(eqn,lastone);
  79.     strcpy(lastone,eqn);
  80.  
  81.     stripspace(eqn);
  82.     strcpy(buff,myfile);
  83.     s1 = strchr(buff,'.');
  84.     if (s1!=NULL) *s1 = 0;
  85.     strcat(buff,".gle");
  86.  
  87.     printf("Output file name to write gle file [%s] ?",buff); gets(outfile);
  88.     if (strlen(outfile)==0) strcpy(outfile,buff);
  89.     printf("Precision of fit required, [1e-4] ?"); gets(buff);
  90.     anstol = atof(buff);
  91.     if (anstol==0) anstol = 1e-4;
  92.     strcpy(equation,eqn);
  93.     lastone_save(eqn,lastdat);
  94.     v_polish(eqn);
  95.     dofit();
  96. }
  97. /*---------------------------------------------------------------------------*/
  98. int xyz_alloc;
  99. double *pntx;
  100. double *pnty;
  101. FILE *fp;
  102. int nxy;
  103. long npnts;
  104. pnt_alloc(int size)
  105. {
  106.     void *d,*dd;
  107.     if (size+10<xyz_alloc) return;
  108.     size = size*2;
  109.     d = farmalloc(size*sizeof(double));
  110.     dd = farmalloc(size*sizeof(double));
  111.     if (dd==NULL) printf("Unable to allocate storage for data %d\n",size*sizeof(double));
  112.     if (d==NULL) printf("Unable to allocate storage for data %d\n",size*sizeof(double));
  113.     if (d==NULL || dd==NULL) exit(1);
  114.     if (xyz_alloc>0) memcpy(d,pntx,xyz_alloc*sizeof(double));
  115.     if (xyz_alloc>0) memcpy(dd,pnty,xyz_alloc*sizeof(double));
  116.     farfree(pntx); farfree(pnty);
  117.     pnty = dd;
  118.     pntx = d;
  119.     xyz_alloc = size;
  120. }
  121. static char buff[2001];
  122. load_data(char *fname,int c1, int c2)
  123. {
  124.     double v;
  125.     char *s;
  126.     FILE *df;
  127.     int i,nd,nc,gotx,goty;
  128.     int done_err;
  129.     done_err = false;
  130.     xmin = 10e10;
  131.     xmax = -10e10;
  132.     pnt_alloc(30);
  133.  
  134.     df = fopen(fname,"r");
  135.     if (df==NULL) {
  136.         printf("Unable to open data file %s \n",fname);
  137.         exit(1);
  138.     }
  139.     nd = 0;
  140.     for (;!feof(df);) {
  141.       if (fgets(buff,2000,df)!=NULL) {
  142.         if (strchr(buff,'!')!=NULL) *strchr(buff,'!') = 0;
  143.         gotx = goty = false;
  144.         s = strtok(buff," \t\n,");
  145.         for (nc=1;s!=NULL;nc++) {
  146.             v = atof(s);
  147.             pnt_alloc(nd);
  148.             if (isdigit(*s) || *s=='-' || *s=='+' || *s=='.') {
  149.             if (nc==c1) {pntx[nd] = v; gotx = true;}
  150.             if (nc==c2) {pnty[nd] = v; goty = true;}
  151.             } else if (*s != '*') gprint("Not a number {%s} \n",s);
  152.             s = strtok(NULL," \t\n,");
  153.         }
  154.         if (gotx && goty) {
  155.             if (pntx[nd]>xmax) xmax = pntx[nd];
  156.             if (pntx[nd]<xmin) xmin = pntx[nd];
  157.             nd++;
  158.         }
  159.       }
  160.     }
  161.     fclose(df);
  162.     nxy = nd;
  163. }
  164. /*---------------------------------------------------------------------------*/
  165. int v_idx[26], naz;
  166. long cpcode[200];
  167. v_polish(char *exp)
  168. {
  169.     int cp,vtype,vidx,ddok,ssok;
  170.     int plen,i;
  171.     vtype = 1;
  172.  
  173.     for (i=0;i<26;i++) v_idx[i] = -1;
  174.     polish(exp,(char *) cpcode,&plen,&vtype);
  175.     var_find_az(v_idx,&naz);
  176.     naz = 1;
  177.     for (i=0;i<'x'-'a'; i++) {
  178.         if (v_idx[i]>0) naz = i+1;
  179.     }    
  180. }
  181.  
  182. double evalfn(double x,double p[]);
  183. double myfunc(double p[]);
  184. int powell(double p[],double **xi, int n,double ftol, int *iter, 
  185.     double *fret, double (*func)());
  186.  
  187.  
  188. double **matrix(int a, int b, int c, int d);
  189. static double pms[28];    
  190. double rr_fit;
  191. dofit()
  192. {
  193.     FILE *fp;
  194.     static double **xi;
  195.     int i,j;
  196.     int niter=0;
  197.     double fret=0;
  198.  
  199.     xi = matrix(1,naz,1,naz);
  200.     for (i=1;i<=naz;i++) for (j=1;j<=naz;j++) xi[i][j] = 0;
  201.     for (i=1;i<=naz;i++) xi[i][i] = 1;
  202.     for (i=1;i<=naz;i++) {
  203.         printf("Initial value for constant %c [1.0] ? ",i+'a'-1); gets(buff);
  204.         if (strlen(buff)==0) strcpy(buff,"1");
  205.         pms[i] = atof(buff);    
  206.     }
  207.     powell(pms,xi,naz,anstol,&niter,&fret,myfunc);
  208.     for (i=1;i<=naz;i++) printf("%c = %g ",i+'a'-1,pms[i]);
  209.     printf("\n");
  210.  
  211.     build_equation(equation,pms,naz);
  212.     stringlower(equation);
  213.     printf("%d Iterations, sum of squares divided by nvalues,  %g\n",niter,(double) fret);
  214.     test_fit();
  215.     printf("y = %s\n",equation);
  216.  
  217.     fp = fopen(outfile,"w");
  218.     fprintf(fp,"size 24 18\n");
  219.     fprintf(fp,"begin graph\n");
  220.     fprintf(fp,"\tsize 24 18\n");
  221.     fprintf(fp,"\ttitle \"y = %s, r^2 = %4.2f%%\" hei .4\n",equation,rr_fit);
  222.     fprintf(fp,"\tdata %s d1=c%d,c%d \n",inputfile,col1,col2);
  223.     fprintf(fp,"\td1 marker dot \n");
  224.     fprintf(fp,"\tlet d2 = %s from %g to %g step %g \n",equation,xmin,xmax,(xmax-xmin)/100.0);
  225.     fprintf(fp,"\td2 line \n");
  226.     fprintf(fp,"end graph\n");
  227.     fclose(fp);
  228. }
  229. test_fit()
  230. {
  231.     int i,j;
  232.     static int count;
  233.     double v,sumy,sumy2,sumfy,sumfy2,meany,meanfy,rr,yy,fyy,sum1,sum2
  234.         ,sum3,y,fy;
  235.  
  236.     sumy = sumy2 = sumfy = sumfy2 = 0;
  237.     for (i=0;i<nxy;i++) {
  238.         var_set(v_idx['x'-'a'],pntx[i]);
  239.         v_eval(&v);
  240.         sumfy = sumfy + v;
  241.         sumy = sumy + pnty[i];
  242.     }
  243.     if (nxy==0) printf("No data points \n");
  244.     meany = sumy/nxy;
  245.     meanfy = sumfy/nxy;
  246.     
  247.     sum1 = 0; sum2 = 0; sum3 = 0;
  248.     for (i=0;i<nxy;i++) {
  249.         var_set(v_idx['x'-'a'],pntx[i]);
  250.         v_eval(&fy);
  251.         y = pnty[i];
  252.         fyy = v-meanfy;
  253.         yy = pnty[i]-meany;
  254.         sum1 = sum1 + (y-fy)*(y-fy);
  255.         sum2 = sum2 + (y-meany)*(y-meany);
  256.     }
  257.     
  258.     rr = 1- (sum1/sum2);
  259.     rr = 100*rr;
  260.     printf("r^2 = %4.2f%%\n",rr);
  261.     rr_fit = rr;
  262.     /* r^2=100*( 1 - (sum of (y-yf)^2)/(sum of (y-ymean)^2) ) */
  263. }
  264. double myfunc(double p[])
  265. {
  266.     int i,j;
  267.     static int count;
  268.     double tot=0,v;
  269.     for (j=0;j<naz;j++) {
  270.         if (v_idx[j]>=0) var_set(v_idx[j],p[j+1]);
  271.     }
  272.     for (i=0;i<nxy;i++) {
  273.         var_set(v_idx['x'-'a'],pntx[i]);
  274.         v_eval(&v);
  275.         tot = tot + pow(fabs(pnty[i] - v),2);
  276.     }
  277.     if (nxy==0) printf("No data points \n");
  278.     tot = tot/nxy;
  279.     if (count/20==count/20.0) {
  280.         printf("%d evaluations, ",count);
  281.         for (i=1;i<=naz;i++) {
  282.             printf("%g ",pms[i]);
  283.         }
  284.         printf(", fit = %g \n",tot);
  285.     }        
  286.     count++;
  287.     /* printf("Evaluate fit = %g\n",tot); */
  288.     return tot;
  289. }
  290. v_eval(double *v)
  291. {
  292.     int cp,vtype,vidx,ddok,ssok,j;
  293.     char outstr[30];
  294.     vtype = 1;
  295.     cp = 0;
  296.     eval(cpcode,&cp,v,outstr,&vtype);
  297. }
  298. gle_abort(char *s)
  299. {
  300.     printf("Abort {%s} \n",s);
  301.     exit(1);
  302. }
  303. getch(){}
  304. stripspace(char *ss)
  305. {
  306.     char *s=ss;
  307.     for (;*s != 0;s++) {
  308.         if (*s != ' ') *ss++ = *s;    
  309.     }
  310.     *ss++ = 0;
  311. }
  312. build_equation(char *expr, double pms[], int naz)
  313. {
  314.     static char inbuff[200];
  315.     static char *tk[500];
  316.     static char tkbuff[500],answer[300];
  317.     static int ntk,ct,i,c;
  318.     char *space_str=" ";
  319.     static char buff[50];
  320.     if (tk[400]==NULL) for (i=0;i<500;i++) tk[i] = space_str;
  321.     answer[0] = 0;
  322.     ntk = 0; 
  323.     token_norm();
  324.     token(expr,tk,&ntk,tkbuff);
  325.     token_space();
  326.     for (ct=1; ct<=ntk; ct++) {
  327.         if (strlen(tk[ct])==1)  {
  328.             c = *tk[ct] - 'A'; 
  329.             if (c>=0 && c<23) {
  330.                 sprintf(buff,"%g",pms[c+1]);
  331.                 strcat(answer,buff);
  332.                 continue;
  333.             }
  334.         }
  335.         if (*tk[ct] != ' ') strcat(answer,tk[ct]);
  336.     }
  337.     strcpy(expr,answer);
  338. }
  339. stringlower(char *s)
  340. {
  341.     for (;*s != 0;s++) *s = tolower(*s);
  342. }
  343.  
  344. /* dummy routines for eval and polish */
  345. g_get_type(){}
  346. g_get_xy(){}
  347. g_measure(){}
  348. tex_xend(){}
  349. graph_xgraph(){}
  350. tex_yend(){}
  351. graph_ygraph(){}
  352. pass_marker(){}
  353. pass_font(){}
  354. pass_color(){}
  355.  
  356. lastone_save(char *s,char *sdat)
  357. {
  358.     FILE *f;
  359.     unlink("fitls.eqn");
  360.     f = fopen("fitls.eqn","w");
  361.     if (f==NULL) {printf("Couldn't save eqn in FITLS.EQN\n"); return;}
  362.     fprintf(f,"%s\n",s);
  363.     fprintf(f,"%s\n",sdat);
  364.     fclose(f);
  365. }
  366. lastone_get(char *s, char *sdat)
  367. {
  368.     FILE *f;
  369.     char *ss;
  370.     f = fopen("fitls.eqn","r");
  371.     if (f==NULL) goto dflt;
  372.     fgets(s,200,f);
  373.     ss = strchr(s,'\n');
  374.     if (ss!=NULL) *ss = 0;
  375.     if (fgets(sdat,200,f)==NULL) goto dflt2;
  376.     ss = strchr(sdat,'\n');
  377.     if (ss!=NULL) *ss = 0;
  378.     fclose(f);
  379.     if (strlen(sdat)==0) goto dflt2;
  380.     return;
  381. dflt:
  382.     strcpy(s,"a+b*x"); 
  383. dflt2:
  384.     strcpy(sdat,"test.dat,1,2");
  385.     return;
  386. }
  387.